#ifndef VARIANCESTRUCTURE_HPP
#define VARIANCESTRUCTURE_HPP
#include<vector>

namespace TermStructure{
    //! A variance structure. Please refer to the doc for more details on the theories
	template<typename T,typename TSize=long>
	class VarianceStructure
	{
        //! Underlying Time Sequence
		TimeSequence<T,TSize> *_times;
        //! utility function
		T Power(T a, T b) const
		{
			return std::pow(a,b);
		}
        //! utility function
		T _integrate(T x0, T x1, T ti, T tj) const
		{
			static T E = M_E;
			if(b==0)
			{
				return (-(x0*(6*(c + d + a*ti)*(c + d + a*tj) - 3*a*(2*c + 2*d + a*(ti + tj))*x0 + 2*Power(a,2)*Power(x0,2))) + x1*(6*(c + d + a*ti)*(c + d + a*tj) - 3*a*(2*c + 2*d + a*(ti + tj))*x1 + 2*Power(a,2)*Power(x1,2)))/6;
			}
			return (Power(a,2)*(-(Power(E,2*b*x0)*(1 + b*(ti + tj - 2*x0) + 2*Power(b,2)*(ti - x0)*(tj - x0))) + 
				Power(E,2*b*x1)*(1 + b*(ti + tj - 2*x1) + 2*Power(b,2)*(ti - x1)*(tj - x1))) - 
				2*Power(b,2)*(2*c*d*(Power(E,b*ti) + Power(E,b*tj))*(Power(E,b*x0) - Power(E,b*x1)) + 
				Power(d,2)*(Power(E,2*b*x0) - Power(E,2*b*x1)) + 2*b*Power(c,2)*Power(E,b*(ti + tj))*(x0 - x1)) + 
				2*a*b*(d*(Power(E,2*b*x0)*(-1 - b*(ti + tj - 2*x0)) + Power(E,2*b*x1)*(1 + b*(ti + tj - 2*x1))) + 
				2*c*(Power(E,b*(tj + x0))*(-1 - b*ti + b*x0) + Power(E,b*(ti + x0))*(-1 - b*tj + b*x0) + 
				Power(E,b*(tj + x1))*(1 + b*ti - b*x1) + Power(E,b*(ti + x1))*(1 + b*tj - b*x1))))/
				(4.*Power(b,3)*Power(E,b*(ti + tj)));
		}
	public:
        //! The Phis
		std::vector<T> Phis;
        //! Other parameters other than Phis
		T a,b,c,d,beta1, beta2;
		//! Set time sequence
		void SetTimeSequence(TimeSequence<T,TSize> &ts)
		{
			_times=&ts;
			Phis.resize(_times->GetSize(), 1);

		}
        //! get volatility for a forward rate i at time t
        /*!
         \param i the index of the forward rate
         \param t the time of the desired volatility
         */
		T GetVolatility(TSize i,T t) const
		{
			if(t>=_times->GetTime(i-1))
				return 0;
			else
				return Phis[i]*((a*(_times->GetTime(i-1)-t)+d)*std::exp(-b*(_times->GetTime(i-1)-t))+c);
		}
        //! get the correlation between two forward rates
		T GetCorrelation(TSize i,TSize j,T t) const
		{
			if(i == j) return 1;
			return std::exp(-beta2*std::abs(_times->GetTime(i-1)-_times->GetTime(j-1)))*(1-beta1)+beta1;
		}
        //! get the number of forward rates
		TSize GetSize() const
		{
			return _times->GetSize();
		}
        //! Integrate the volatility from time fromT to ToT.
		T VolatilityIntegrator(T fromT, T toT, TSize i, TSize j) const
		{
			if(toT >= this->_times->GetTime(i-1)) toT = this->_times->GetTime(i-1);
			if(toT >= this->_times->GetTime(j-1)) toT = this->_times->GetTime(j-1);
			if(toT <= fromT) return 0;
			else return this->_integrate(fromT, toT, this->_times->GetTime(i-1), this->_times->GetTime(j-1)) * this->GetCorrelation(i,j,0) * this->Phis[i] *this->Phis[j];
		}
	};

	//! Initiate an instantaneous covariance matrix from a variance structure
	template <typename T, typename TVarianceStructure, typename TSize = long>
	class VarianceStructureMatrix     
	{
	protected:
		TVarianceStructure *_underlyingStructure;
	public:
		T CurrentTime;
		void SetVariaceStructure(TVarianceStructure &underlying)
		{
			_underlyingStructure=&underlying;
		}


		T operator()(TSize i, TSize j) const
		{
			return std::sqrt(_underlyingStructure->GetVariance(i+2, CurrentTime) * _underlyingStructure->GetVariance(j+2, CurrentTime)) * _underlyingStructure->GetCorrelation(i+2,j+2,CurrentTime);
		}

		TSize GetLength(int i) const
		{
			return _underlyingStructure->GetSize()-2;
		}
	};
}


#endif